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Abstract 

The distributed-relaxation multigrid and defect- 
correction methods are applied to the two- 
dimensional compressible Navier-Stokes equations. 
The formulation is intended for high Reynolds num- 
ber applications and several applications are made at 
a laminar Reynolds number of 10,000. A staggered- 
grid arrangement of variables is used; the coupled 
pressure and internal energy equations are solved 
together with multigrid, requiring a block 2x2 ma- 
trix solution. Textbook multigrid efficiencies are at- 
tained for incompressible and slightly compressible 
simulations of the boundary layer on a flat plate. 
Textbook efficiencies are obtained for compressible 
simulations up to Mach numbers of 0.7 for a viscous 
wake simulation. 

Introduction 

Computational fluid dynamics (CFD) is becom- 
ing a more important part of the complete aircraft 
design cycle because of the availability of faster com- 
puters with more memory and improved numerical 
algorithms. The cruise shapes of transport aircraft 
are designed to minimize viscous and shock wave 
losses at transonic speeds and computational meth- 
ods for these flows are reasonably well in hand. Sim- 
ulations of off-design performance associated with 
maximum lift, buffet, and flutter and the deter- 
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min at ion of stability and control derivatives, in- 
volving unsteady separated and vortical flows with 
stronger shock waves, demand significant comput- 
ing resources since Reynolds-averaged Navier-Stokes 
(RANS) methods are needed. The turnaround time 
of these computations for high-Reynolds flows over 
complex geometries is too long to impact the design 
cycle and the turbulence models for separated flows 
have a high degree of variability. 

The current RANS solvers with multigrid require 
on the order of 1500 residual evaluations to con- 
verge the lift and drag to one percent of their fi- 
nal values for wing-body geometries near transonic 
cruise conditions. Complex geometry and complex 
physics simulations require many more residual eval- 
uations to converge, if indeed convergence can even 
be attained. It is well-known for elliptic problems 
that solutions can be attained optimally using a full 
multigrid (FMG) process in far fewer, on the order of 
2-4, residual evaluations. An optimally convergent 
method is defined by Brandt 1 ' 2 as textbook multi- 
grid efficiency (TME), meaning the solutions to the 
governing system of equations are attained in a com- 
putational work which is a small (less than 10) mul- 
tiple of the operation count in the discretized system 
of equations. Thus, there is a potential gain of two 
orders of magnitude in operation count reduction if 
TME could be attained for the RANS equation sets. 
This possible two order of magnitude improvement 
in convergence represents an algorithmic floor since 
it is unlikely that faster convergence for these non- 
linear equations could be attained. This algorithmic 
acceleration, however, coupled with further increases 
in computational speed can open up avenues and ac- 
celerate progress in many areas, including: the ap- 
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plication of steady and time-dependent simulations 
in the high-lift, off-design, and stability and control 
areas; the usage of RANS solvers in the aerodynamic 
and multidisciplinary design areas; and the develop- 
ment of improved turbulence models. 

The RANS equation sets are a system of coupled 
nonlinear equations which are not, even for subsonic 
Mach numbers, fully elliptic, but contain hyperbolic 
partitions. Brandt 1 has summarized the progress 
and remaining barriers in TME for the equations 
of fluid dynamics. Although TME demonstrations, 
especially for incompressible simulations, have been 
completed, the distributed relaxation approach has 
not been used widely for large-scale simulations. 
Current methodologies use a block-matrix relaxation 
and/or a time-dependent approach to solve the equa- 
tions; significant improvements have been demon- 
strated using multigrid approaches, but the meth- 
ods are not optimally convergent. An example is 
the plane solver of Thomas, et al. 3 which gives fast 
multigrid performance, but at a cost of block matrix 
solutions of size m, where m is the number of conser- 
vation equations at a grid point. A subtle disadvan- 
tage of this block-matrix approach is that classical 
Gauss-Seidel relaxation methods, well suited for the 
constituent hyperbolic and elliptic partitions in the 
equations, cannot be applied stably in subsonic flow. 

The distributed relaxation approach of Brandt de- 
composes the system of equations into separable, 
many times scalar, pieces that can be treated with 
optimal methods. The purpose of this paper is to 
apply distributed relaxation in a defect correction 
setting to the compressible Navier-Stokes equations. 
A staggered- grid discretization is used which sim- 
plifies the elements of the calculation, although this 
is not an essential ingredient. The method follows 
closely the original derivation of Brandt, 1,2 dating 
to 1984. 


Distributed Relaxation 

The equations for the time-dependent conserva- 
tion of mass, momentum, and energy can be written 
as 


<9 t Q + R = 0, (1) 

where Q = (p, pu, pv, e) T and R(Q) is the spatial 
divergence of a vector function representing convec- 
tion and viscous/heat transfer effects. The quantity 
e is the total energy per unit volume, € is the internal 
energy per unit mass, and the perfect gas equation 
of state relating pressure, p, to density, p, and e can 
be written as 


P = ( 7 - l)[e - p{u 2 + u 2 )/2] = (7 - 1 )pc, (2) 

where j is the ratio of specific heats. A set of non- 
conservative equations can be cast in terms of prim- 
itive variables q = (u,v,p, e) T by multiplying the 
equations by the Jacobian matrix given as 


dq_ 

dQ 


= S 


-1 1 0 

-1 0 1 

(u 2 -j- v 2 )/ 6 ! — u —v 

— e -f (u 2 + t> 2 )/2 — u —v 


0 

0 

1 

1 


( 3 ) 


S EE 


0 

p 

0 

0 


0 

0 

7-1 

0 


0 

0 

0 

p 


( 4 ) 


The primitive equations are then 


ft,+ ||R = 0. 


( 5 ) 


These equations may be linearized about a given 
state variable and written in delta form, £q = 


1^+1 


as 


= < 6 > 

where in the matrix operator L, only the principal 
terms at the viscous and inviscid scales are retained. 
Both scales are essential for high Reynolds number 
simulations: the inviscid scales over most of the flow- 
field and the viscous scales in the thin viscous layers 
near bodies and in their wakes. Since we are in- 
terested here in solutions to the steady-state equa- 
tions, we drop the time derivative. Additionally, the 
thin-layer approximation, in which only the viscous 
terms associated with variations in the coordinate 
normal to the body are retained, can further sim- 
plify the terms in L. This approximation is widely- 
used for high Reynolds number applications and is 
used for the numerical calculations below, although 
the development does not rely on this approxima- 
tion. Also, the formulation allows a conservative 
discretization of R although, for the results below, 
a nonconservative form for R is used. 

The distributed relaxation method replaces <Jq by 
Mtfw so that the resulting matrix LM becomes a 
diagonal or lower triangular matrix, as 
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LMSw = ~^|REE-r. (7) 

The diagonal elements of LM are composed ideally 
of the separable components of the determinant of 
the matrix L and represent the elliptic or hyperbolic 
features of the equation. 

The approach yields fast convergence for both 
steady and unsteady simulations if the constituent 
scalar diagonal operators in LM are solved with fast 
methods, such as multigrid for the elliptic partitions 
and upwind methods for convection. Brandt 1 has 
derived a set of matrices M for the equations of fluid 
dynamics which provide a convenient lower triangu- 
lar form based on a discretization of the equation of 
state. A variation of that formulation, in which the 
equation of state is an algebraic equation, is used 
below. The Appendix gives the matrices for the in- 
compressible equations. 

Defect Correction 

Since Eq. (7) is written in delta form, it is natural 
to consider defect correction for the update, namely 
a lower-order discretization of the left side of the 
equations in order to simplify the construction of 
the implicit equations. The most common approx- 
imation for the driver operator is a first-order dis- 
cretization for the convective and pressure (inviscid) 
contributions. For this approach, we can write the 
implicit scheme as 


[LM] d (Sw = —iv (8) 

where the subscripts t and d denote some desired 
“target” and “driver” schemes on the right and left 
sides, respectively, of the equation. 

The convergence of defect correction is known to 
be fast for elliptic equations. For hyperbolic equa- 
tions, the convergence may be slow, even if the im- 
plicit equations are solved exactly. This slowdown 
was pointed out by Brandt 5 and discussed in some 
detail in Brandt and Yavneh, 4 Thomas, et. al, 3 and 
Diskin and Thomas/ For the second-order upwind- 
biased discretization corresponding to k = 0, defined 
subsequently, the asymptotic convergence rate is ap- 
proximately 0.5 per defect-correction iteration. For 
a first-order driver approximation, the initial conver- 
gence may be quite slow; the number of iterations 
to get into the asymptotic convergence regime might 
grow on fine grids and for certain frequencies at the 
rate h~ U 3 . This worst case behavior is usually not 


that harmful for steady-state calculations with mod- 
erate frequency content of the information imposed 
at inflow. 

With a second-order accurate behavior of the 
driver operator such as the fully-upwind second- 
order operator (k = — 1), the convergence is quite 
fast, within three iterations at most, although on 
very coarse meshes, in which neither the driver nor 
the target operators retain discretely an order prop- 
erty over the whole of the domain, some large growth 
of errors may occur. An efficient half-space Fourier 
mode analysis is given in Diskin and Thomas 7 which 
predicts quite closely the actual solution behavior for 
scalar equations. 

Inviscid Compressible Equations 

The inviscid compressible equations take a rela- 
tively simple form with the choice of primitive vari- 
ables considered here. Using a nonconservative dis- 
cretization of the steady equations, the update equa- 
tion is defined as 


LSq — r — -(f -f Lq), (9) 


where f = 0, the matrix L is 


Q 0 jd x 0 ' 

0 Q jdy 0 

pc 2 d x pc 2 d y Q 0 1 

c ~d x c ydy 0 Q 


( 10 ) 


and Q = ud x + vd y . The last equation is decoupled 
from the first three, representing the convection of 
entropy along a streamline. The determinant of the 
matrix of operators, 


g 2 [Q 2 -c 2 A], 

corresponds to an elliptic partition, represented by 
the full potential equation, and two convective par- 
titions, generally recognized as the convection of en- 
tropy and vorticity. Taking the distribution matrix 
M as the cofactors of the third row of L divided by 
their common factor, 


M = 


1 0 -ifl* 0 

0 1 ~\dy 0 

0 0 Q 0 

0 0 0 1 


( 11 ) 


the equations to solve for the “ghost” Sw variables 
are 
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LM<Jw = -r, (12) 


or 


QSw i 
QSw 2 


= -n 


{—^Q 2 + \) 6 w 3 = 
^QSw 4 = 

c 2 


~r 2 

4*^2 r 3 + p(d x Swi + d y Sw 2 ) 

— ~~2 r 4 ~ d x 8w i + d y Sw 2 ) 

+ — A8w 3 
PI 


The last equation could be replaced with 


QSq^ — — r 4 - -~(d x 8u + d y 8v) 
for inviscid flow. Thus the update to q is 


!l n+1 =U n +6 qi 


V n + l Sg2 

p n+l =^+Sq 3 

e n +l = 6 n + Sg4 


u n -f 8wi d x Sw 3 

P 

v n -f 8w 2 — -d y Sw 3 

p n + Q8w 3 
e n -f 8W4 


The staggered- grid discretization used here is usual: 
p,e, p defined at the interiors of the grid, u defined at 
the cell interfaces tangent to the y— or k— direction, 
and v defined at the cell interfaces tangent to the 
or j— direction. The discrete scheme with such a 
staggered-grid arrangement of variables can be de- 
scribed as 


Q h o 

_ 0 Q* # 

pc 2 d h x pc 2 d£ Q% 

. T 9 " c i d y 0 

where d x and d x are 2-point centered differences on 
the staggered grid. The operator Q h is an upwind- 
biased convection operator, defined in terms of 
translation operators 7) and 7 * ■ (TjU j k = Uj + \ k )- 
as 


0 

0 

0 

Q h 


(13) 


Q h = 

D-(Tj) = 

D + (Tj) = 

= 

The coefficients for k G [— 1, 1] are 

{ci,c 2 ,c 3 ,c 4 } = ^{1 — Kj 3k — 5, 3(1 - /c), 1 + k} 
and for first-order accuracy are 


tr „ x u 

D-(Tj) + 


(M 


[hx)j 


D + (Tj) 

-D + (Tk) 


?)' V 

+ /,'\ D ~(Tk) + 7T T 

(hy)k (hy), 

ciTj 2 + c 2 7^ ^ + C 3 + c^Tp 
- Cl T/ 2 - c 2 Tp - c 3 - C4T- 1 

(«± |«D/2 


{ci,c 2 ,c 3 ,c 4 } = {0, -1, 1, 0}. 

The convection operator in the pressure equation is 
differenced with a downwind operator, i.e., 


<58 = 


W D+ra W‘ (Tf) 

+7^0 + ™+7 frD-(T.) 


so that the determinant of the discrete matrix L h , 


[Q h ?[Q h Q h o - c 2 A h ], 

has a smaller and more symmetrical stencil than 
if the upwind operator were used everywhere. For 
grid-aligned flow, this full potential operator is a cen- 
tered 5-point discretization of the Prandtl-Glauert 
equation. With this formulation, the low-Mach 
number limit is not a problem, as the usual five-point 
Laplacian is obtained for the w 3 variable, analogous 
to the purely incompressible form. 

Compressible Navier-Stokes Equations 

The contributions to the principal terms of the 
matrix of operators for the viscous and heat con- 
duction terms arise from the momentum and energy 
contributions as in Eq. (3), since there is no influ- 
ence of viscosity in the continuity equation. The 
two major complications are that the heat conduc- 
tion term couples into both the pressure and inter- 
nal energy equations and the viscous terms in the 
momentum equations involve cross- derivative terms. 
With the nonconservative discretization considered 
here assuming constant viscosity and heat conduc- 
tion coefficients, the update equation can be written 
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as Eq. (9) with f = -(0, 0, (7 - 1)4>, p$) T where the 
viscous dissipation term $ is 


< 2 > = p[2(d x u) 2 -\-2(dyv)’ i - l r(d x v- ! rdyu)‘ 2 }->r\[{d x u+dyV ) 2 ] 


the matrix L is 


L = 


Qv - )dx 

pc 2 d x 

c -d 

^ u x 


p°xy 

Qu — ^dyy 

pC 2 0y 

C -d 

1 °y 


jdy 

Q 

0 


0 

0 

—(7 — l)kA 

Qkj p 

(14) 


and Q x = Q — yA. Nondimensionalizing by density 
and speed of sound and using the Stokes hypothesis 
for the bulk viscosity term, the coefficients become 
v = p/p — Moo / (Rep), k = Mooi/(RePr) and A = 
A -b p = /i/3. 

The distribution matrix M is developed similarly 
to that for the inviscid equations with some addi- 
tions that eliminate the cross-derivative terms, as 
below 


M = 


1 0 -ia, 0 

1 ~ l p dy 0 

Xd x Xd y Q x 0 

0 0 0 " 1 


(15) 


The matrix LM now has only diagonal entries of Q v 
in the first two rows. The last two equations are still 
coupled, requiring a block 2x2 solution as below, 



Sw 3 

Qk/p 

Sw 4 


where 



(16) 


9 3 


9 4 


^2 r 3 + p(d x Sw 1 + d y Sw 2 ) 
+~^Q(d x Swi + dySu>2 ) 


1 i 

2 ^*4 (d x Swi + d y Sw 2 ) 

c 


In the high Reynolds number simulation considered 
here, the term (—^QQ V+ % + A) in Eq. (16) is ap- 

proximated as the full potential operator associated 


with the inviscid equations. The discrete approxi- 
mation of the terms above are straightforward ad- 
ditions to the inviscid compressible terms, since the 
viscous terms are differenced centrally. 

In the calculations below, the thin-layer approxi- 
mation is made for the viscous terms, i.e., only the 
derivative terms in the p-direction are retained. The 
Swi and Sw 2 equations (the momentum equations at 
constant pressure) are solved with a tridiagonal y- 
line marching algorithm. The equations for Sw 3 and 
Svl >4 are solved with two passes of a V(2,l) multi- 
grid cycle with line Gauss-Seidel iteration using four 
multigrid levels; a correction scheme (CS) multigrid 
was used in which the coefficients of the block ma- 
trix terms were interpolated to cell center locations 
using bilinear interpolation. The coarsest grid level 
was solved exactly. A smoothing was done with a 
point Gauss-Seidel scheme before transferring resid- 
uals to coarser meshes. 

Computational Results 

Inviscid Linearized Flow 

The linearized flow over a bump in a channel 
was computed for a computational domain extend- 
ing from x = 0 to x = 3 and y = 0 to y = 1 on a uni- 
form grid of mesh size h x — h y . A sine-squared pro- 
file extending from x — 1 to x — 2 of height 0.1 that 
of the channel height was imposed with linearized 
boundary conditions, v/uoo = dy/dx at y — 0. 

The boundary conditions for the primitive and 
ghost variables require some discussion. The residu- 
als at the horizontally oriented boundaries are taken 
to be zero, corresponding to an imposed boundary 
condition of tangency. In general, satisfaction of 
these residuals to zero is obtained by adjusting the 
pressure across the boundary to satisfy the normal- 
momentum equation. The boundary condition im- 
posed on the Sws ghost variable along the top and 
bottom of the domain is d y Sw 3 = 0, corresponding 
to the prescribed normal velocity condition, Sv = 0. 

The residuals and u at the vertically oriented up- 
stream boundary are unknowns in the solution; the 
variables p, e and u at points just upstream of the 
boundary are taken as the imposed boundary con- 
ditions; the upstream boundary conditions imposed 
in the ghost- variable solutions are as follows: 


Swi = 0 ; Sws =3 0 ; Sw 4 = 0 . 


corresponding to imposed values of u, p, e. The value 
of v upstream of the boundary can be determined by 
imposing Sw 2 = 0 and then correcting it at values 
outside the domain using the update formulas. 
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5r 

II 

3- 



c P , 

M co = 0 

Cp, 

= 0.01 

c p , 

Moo = 0.5 

1/4 

-0.0058 

-0.0058 

-0.0074 

1/8 

-0.1767 

-0.1767 

-0.2116 

1/16 

-0.3164 

-0.3165 

-0.3730 

1/32 1 

-0.3650 

-0.3651 

-0.4274 

1/64 

-0.3785 

-0.3786 

-0.4423 


Table 1. Linearized pressure coefficients at x = 1.5 
for the flow over a sine-squared bump in a channel. 


The pressures at locations just downstream of the 
downstream boundary are specified as freestream; 
all of the other values are obtained through integra- 
tion of the update equations with the upwind-biased 
convection operator. The downstream boundary 
conditions imposed on the ghost variable Sw 3 is 
d x Sw 3 — 0, corresponding to the prescribed pres- 
sure condition, Sp = 0 . Thus, in summary, Dirich- 
let conditions are imposed for Sw at the upstream 
boundary and Neumann conditions for Sw 3 at all 
other boundaries. 

The exact discrete solution for both incompress- 
ible and compressible flow is obtained in one iter- 
ation starting from freestream conditions if (1) we 
solve the Sw 3 equation to machine zero using con- 
ventional multigrid methods and (2) we use first- 
order differencing. We are not recommending that 
residuals be reduced to machine zero; rather this 
approach does provide a convenient check of the 
implementation of the algorithm on the computer, 
since arbitrarily large perturbations in q can be 
damped/con vected from the domain in one iteration. 
The ghost variable Sw% in the domain takes the role 
of the perturbation potential in an irrotational isen- 
tropic flow. The values of u and v are obtained to 
second order. The value of p is determined by the ac- 
curacy of the convection operator. Fast convergence 
and second-order accuracy is obtained if we use de- 
fect correction with the second-order convection op- 
erator (k = 0). The values of C p ~ — 2( « — i*oo) at 
x = 1.5 obtained on a sequence of grids are shown in 
Table 1. Second-order accuracy is obtained and the 
results obey the expected Prandtl-Glauert scaling of 
the pressure coefficient with l/\/(l — M<^); the C p 
decrease for M 0 0 — 0.5 is 1.17, as compared to 1.16 
from the Prandtl-Glauert scaling. 

Viscous Flat Plate 

The viscous flow over a flat plate was computed 
for the same computational domain as above at a 
Re= 10,000 based on the height of the channel. The 
grid was stretched in the y— direction with a stretch- 
ing factor of 1.06 on the finest mesh of N x xN y = 


193x65. No-slip, adiabatic wall conditions are pre- 
scribed from x — 1 to x = 2 along the lower bound- 
ary and symmetry conditions upstream and down- 
stream of those points; a wake profile develops down- 
stream of the trailing edge x — 2. Pressure was pre- 
scribed at the downstream boundary and tangency 
along the upper wall. 

The ghost-variable boundary conditions are the 
same as the channel flow simulation except in the 
region of the flat plate. For the internal energy, 
the gradient normal to the boundary is zero for 
both the adiabatic and symmetry conditions, so that 
d y Sw 4 = 0 is imposed. Since the residuals and v 
are prescribed to be zero at the lower boundary, the 
boundary conditions for the Sw 3 variable is taken 
as dySw 3 = 0. Thus, boundary conditions for the 
ghost variables in this simulation are the same as 
the channel flow simulation. The primitive variable 
velocities outside the lower boundary in the region of 
the plate are chosen to satisfy the no-slip condition 
at the plate, i.e. u(— y) — ~u(y);v(—y) = — v(y ). 

Convergence could be attained in this simulation 
for incompressible and slightly compressible flow. 
An FMG cycle was used. The value of Cf at x = 1.5, 
midway down the plate, is a sensitive measure of con- 
vergence. First-order solutions were obtained which 
decreased the L^-norm of the maximum residual 4-5 
orders of magnitude on each of the three finer meshes 
over 20 iterations, corresponding to a convergence 
rate of roughly 0.6 per fine-grid iteration for both in- 
compressible flow and Moo =0.1. The skin friction 
values converged to within discretization error in 2-4 
cycles. For second-order accuracy of the convection 
operator, k = 0, skin friction results are tabulated 
in Table 2 for the series of meshes and for 2,4, 10, and 
20 fine-grid iterations. The convergence is quite fast; 
solutions are obtained to within truncation-error ac- 
curacy within a few multigrid cycles. The f/2-norm 
of the maximum residual was reduced 3 orders of 
magnitude on each of the three finer meshes in 20 
iterations, corresponding to a convergence rate of 
roughly 0.7 per fine-grid iter at ion, again for both in- 
compressible flow and Moo = 0.1. 

At higher Mach numbers, convergence could not 
be attained. The residuals remained quite large 
and the solution eventually diverged near the lead- 
ing edge of the plate. This behavior is not entirely 
unexpected since we expect the distributed relax- 
ation approach to be augmented with procedures to 
locally reduce the residuals near boundaries, before 
relaxing the interior equations. This local relaxation 
near boundaries is required to guarantee convergence 
within a few cycles even for elliptic equations. 6 
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N x xNy 

Iterations 

Cj, 

Moo = 0 

Cf, 

Moo = 0.1 

25 x 9 

2 

0.973E-02 

0.111E-01 

25 x 9 

4 

0.392E-02 

0.685E-02 

25 x 9 

10 

0.490E-02 

0.574E-02 

25 x 9 

20 

0.482E-02 

0.566E-02 

49 x 17 

2 

0.694E-02 

0.756E-02 

49 x 17 

4 

0.780E-02 

0.831E-02 

49 x 17 

10 

0.813E-02 

0.875E-02 

49 x 17 

20 

0.814E-02 

0.877E-02 

97 x 33 

2 

0.924E-02 

0.934E-02 

97 x 33 

4 

0.937E-02 

0.949E-02 

97 x 33 

10 

0.931E-02 

0.948E-02 

97x33 

20 

0.931E-02 

0.948E-02 

193 x 65 

2 

0.970E-02 

0.970E-02 

193 x 65 

4 

0.967E-02 

0.971E-02 

193 x 65 

10 

0.960E-02 

0.967E-02 

193 x 65 

20 

0.959E-02 

0.967E-02 


Table 2. Convergence of the skin friction value 
halfway down the viscous flat plate for incompress- 
ible and slightly compressible flow; Re — 10,000. 


Viscous Wake 

In order to investigate the behavior of the viscous 
algorithm for higher Mach numbers in a smoother 
flow, the viscous flow in a developing wake profile 
was considered. The simulation was similar to that 
above except that a wake deficit with uniform pres- 
sure and internal energy was applied at the inflow 
boundary. The inlet profile at x — —h x was pre- 
scribed according to the exact incompressible far 
wake solution, as 

u/uco = 1 - w d {x + l)~ 1/2 exp( 2 ) 

4 (x -f 1 y 

where wj = 0.5. The minimum velocity was moni- 
tored at x — 1.5, a location midway in the domain, 
as a measure of convergence. The boundary con- 
ditions for the ghost variables are the same as the 
channel flow simulation. 

Convergence could be attained from incompress- 
ible flow up to the highest Mach number investi- 
gated, Moo = 0.7. First-order solutions were ob- 
tained which decreased the i^-norm of the maxi- 
mum residual 5 orders of magnitude on each of the 
three finer meshes over 10 iterations, corresponding 
to a convergence rate of roughly 0.3 per fine-grid it- 
eration at all Mach numbers The wake deficit values 
converged to within discretization error in 2 itera- 
tions. For second-order accuracy of the convection 


N%xNy 

Iterations 

'Umin 7 
M eo = 0 

l^min j 

Moo = 0.7 

25 x 9 

1 

0.7149 

0.7149 

25 x 9 

2 

0.7116 

0.7141 

25 x 9 

4 

0.7105 

0.7112 

25 x 9 

10 

0.7105 

0.7112 

49 x 17 

1 

0.7474 

0.7478 

49 x 17 

2 

0.7483 

0.7487 

49 x 17 

4 

0.7482 

0.7485 

49 x 17 

10 

0.7481 

0.7485 

97 x33 

1 

0.7520 

0.7523 

97 x 33 

2 

0.7531 

0.7534 

97 x 33 

4 

0.7530 

0.7533 

97 x 33 

10 “ 

0.7530 

0.7533 

193 x 65 

1 

0.7529 

0.7532 

193 x 65 

2 

0.7532 

0.7535 

193 x 65 

4 

0.7531 

0.7535 

193 x 65 

10 

0.7531 

0.7534 


Table 3. Convergence of the minimum velocity at 
x = 1.5 for incompressible and compressible wake 
flows; Re = 10, 000. 


operator, k — 0, minimum velocity values are tab- 
ulated in Table 3 for the series of meshes used in 
the FMG cycle and for 1,2,4, and 10 fine-grid itera- 
tions. The solution converges to within truncation- 
error accuracy within 2 multigrid cycles. The L^~ 
norm of the maximum residual was reduced 3 or- 
ders of magnitude on each of the three finer meshes 
in 10 iterations, corresponding to a convergence rate 
of roughly 0.5 per fine-grid iter at ion, again for both 
incompressible flow and Moo = 0.7. 

Concluding Remarks 

The distributed-relaxation multigrid and defect- 
correction methods are applied to the compress- 
ible Navier-Stokes equations. The formulation is in- 
tended for high Reynolds number applications and 
several applications are made at a laminar Reynolds 
number of 10,000. Although deemed not essential 
to the methodology, both a staggered-grid arrange- 
ment of variables and nonconservative forms for the 
governing equations have been used. 

The compressible flow algorithm solves a multi- 
grid problem for the pressure and internal energy 
equations which necessitates a local block 2x2 ma- 
trix solution at every grid point. Since we obtain 
solutions to within truncation error in 2-4 fine-grid 
iterations and the additional operations to solve for 
the ghost variables are on the order of an additional 
residual evaluation, textbook multigrid efficiencies 
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are attained for incompressible and slightly com- 
pressible simulations of the boundary layer on a flat 
plate. Difficulties were encountered for higher Mach 
numbers; the calculation can and should be aug- 
mented with local procedures to reduce the resid- 
uals near boundaries, since this approach is required 
for general conditions even for the elliptic equations. 
Textbook efficiencies are obtained for compressible 
simulations up to Mach numbers of 0.7 for a viscous 
wake simulation. 
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L == 


Qv 0 d x 0 
0 Q v d y 0 
d x by 0 0 

0 0 0 Q k 


f = — (0, 0, 0, <F) T , and A = 0. Nondimension- 
alized by density and velocity, the viscosity and 
heat conduction terms are v = p/p = 1/Re and 
k — 7 /(RePr). The distribution matrix M is deter- 
mined by the cofactors of the third row of L divided 
by their common factor, as 


1 0 0 

0 1 -dy 0 

0 0 Q v 0 

0 0 0 1 
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